#!/usr/bin/Rscript

#this program is a testing script for certain correlations between rainfall 
# and inflow.  This script fits a LWP model and resamples the residuals to 
# get an estimate of the confidece intervals.
# This was an ititial attempt and the results from this script are not used
# anywhere to produce results.
#
#Author: Cameron Bracken 
#        cwb12@humboldt.edu
#
#Last Modeified: 7 Dec 2008

source('mylag.R')
source('callCollect.R')
library(locfit)

file = '../overlap/eka-ppt-inch-arc-inflow-mgd-overlap.txt'
baseflow =  0
resample = T
nsims = 100
nn = 10
collect = 3

#Set a lag to shift the dependent variable back forward units in time 
#eg: if y=inflow and x=rainfall setting a lag of 1 will correlate tomorrows
#    inflow with todays rainfall
lag = 0
alpha=.5

    #read in the data 
data = read.table(file,fill=T,header=T)
names = names(data)
data = as.matrix(data)
data[data==-9999]=NA
x = data[,5]
y = data[,4] - baseflow
xlab = names[5]
ylab = names[4]

#Set a logical filter 
filter <- which(x>.3)


time = data[,1]*10000+data[,2]*100+data[,3]
time = as.POSIXlt(as.character(time), tz = "PDT", format = "%Y%m%d")

oo = (1:length(y))[filter]

time = time[filter]
x = x[filter]
y = y[filter]

co = callCollect(x=x, y=y, collect = collect, plot = T)

time = time[collect:length(as.double(time))]
x = co$x
y = co$y

quartz()
plot(x,mylag(y,lag),xlab = xlab, ylab = ylab,cex=.7)

cat('Correlation coefficient:',cor(x,mylag(y,lag),use='complete.obs'),'\n')
fit = lm(mylag(y,lag)~x)
lfit = locfit(mylag(y,lag)~lp(x,alpha=alpha))
abline(fit,col='steelblue')
lines(lfit,col='dodgerblue')

pred = vector('numeric',length(x))

if(resample){ 
    noise = vector('numeric',nsims)
    upper = vector('numeric',length(x))
    lower = vector('numeric',length(x))
}

for(i in 1:length(x)){
    pred[i] = predict(fit,newdata=data.frame(x=x[i]))#
    if(resample){

        d = order(abs(x[i]-x))[2:(nn+1)]
        se = sd(fit$residuals[d])

        for(j in 1:nsims){
            noise[i] = sample(fit$residuals[d],1)
        }
        upper[i] = pred[i] + quantile(noise,.95,na.rm=T)
        lower[i] = pred[i] + quantile(noise,.05,na.rm=T)
    }
}

window = 200:300
quartz(width=14)
plot(y[window],type='l',lwd=.5)
lines(pred[window],col='red',lwd=.5)
lines(upper[window],lty='dashed',lwd=.5)
lines(lower[window],lty='dashed',lwd=.5)



